sink("survey.outcomes-psrm-log.txt",append=F,type="output")

library(plyr)
load('individual.use.election-psrm.RData')

###
#Graphical Comparison
###

keep.ids <- which(!is.na(individual.use.election$candidate.diff.pre) & !is.na(individual.use.election$candidate.diff.post))
keep.frame <- individual.use.election[keep.ids,]

candidate.assessment.bytime <- ddply(keep.frame,.(rep.partisan.bin),summarise,pre.trump=mean(candidate.diff.pre,na.rm=TRUE),post.trump=mean(candidate.diff.post,na.rm=TRUE), pre.trait=mean(candidate.trait.pre, na.rm=TRUE), post.trait=mean(candidate.trait.post, na.rm=TRUE), pre.emotion=mean(candidate.feel.pre, na.rm=TRUE), post.emotion=mean(candidate.feel.post, na.rm=TRUE))

candidate.assessment.bytime.content <- na.omit(ddply(keep.frame,.(rep.partisan.bin.content),summarise,pre.trump=mean(candidate.diff.pre,na.rm=TRUE),post.trump=mean(candidate.diff.post,na.rm=TRUE), pre.trait=mean(candidate.trait.pre, na.rm=TRUE), post.trait=mean(candidate.trait.post, na.rm=TRUE), pre.emotion=mean(candidate.feel.pre, na.rm=TRUE), post.emotion=mean(candidate.feel.post, na.rm=TRUE)))

#Plot Results
pdf(file='figure-4.pdf',height=8,width=8)

par(mar=c(4.1,4.1,1.2,.5),mfrow=c(3,2))
plot(x=1:5,y=candidate.assessment.bytime[,2],type='l',lwd=3,col='black',ylab='Trump - Clinton Feeling Therm',xlab='',las=1,ylim=c(-80,80),main='Audience-Based Exposure Measure',cex.lab=1.2,cex.axis=1.2)
points(x=1:5,y=candidate.assessment.bytime[,2],pch=16,cex=1.5)
lines(x=1:5,y=candidate.assessment.bytime[,3],type='l',lwd=3,col='gray')
points(x=1:5,y=candidate.assessment.bytime[,3],pch=16,col='gray',cex=1.5)
legend('topleft',legend=c('Survey Wave 1','Survey Wave 2'),col=c('black','gray'),lwd=c(4,4), cex=1.2)

plot(x=1:5,y=candidate.assessment.bytime.content[,2],type='l',lwd=3,col='black',ylab='',xlab='',las=1,ylim=c(-80,80),main='Content-Based Exposure Measure',cex.lab=1.2, cex.axis=1.2)
points(x=1:5,y=candidate.assessment.bytime.content[,2],pch=16,cex=1.5)
lines(x=1:5,y=candidate.assessment.bytime.content[,3],type='l',lwd=3,col='gray')
points(x=1:5,y=candidate.assessment.bytime.content[,3],pch=16,col='gray',cex=1.5)

plot(x=1:5,y=candidate.assessment.bytime[,4],type='l',lwd=3,col='black',ylab='Trump - Clinton Trait Rating',xlab='',las=1,ylim=c(-15,15),main='',cex.lab=1.2, cex.axis=1.2)
points(x=1:5,y=candidate.assessment.bytime[,4],pch=16,cex=1.5)
lines(x=1:5,y=candidate.assessment.bytime[,5],type='l',lwd=3,col='gray')
points(x=1:5,y=candidate.assessment.bytime[,5],pch=16,col='gray',cex=1.5)

plot(x=1:5,y=candidate.assessment.bytime.content[,4],type='l',lwd=3,col='black',ylab='',xlab='',las=1,ylim=c(-15,15),main='',cex.lab=1.2, cex.axis=1.2)
points(x=1:5,y=candidate.assessment.bytime.content[,4],pch=16,cex=1.5)
lines(x=1:5,y=candidate.assessment.bytime.content[,5],type='l',lwd=3,col='gray')
points(x=1:5,y=candidate.assessment.bytime.content[,5],pch=16,col='gray',cex=1.5)

plot(x=1:5,y=candidate.assessment.bytime[,6],type='l',lwd=3,col='black',ylab='Trump - Clinton Emotional Reaction',xlab='Republican News Exposure (Low to High)',las=1,ylim=c(-15,15),main='',cex.lab=1.2, cex.axis=1.2)
points(x=1:5,y=candidate.assessment.bytime[,6],pch=16,cex=1.5)
lines(x=1:5,y=candidate.assessment.bytime[,7],type='l',lwd=3,col='gray')
points(x=1:5,y=candidate.assessment.bytime[,7],pch=16,col='gray',cex=1.5)

plot(x=1:5,y=candidate.assessment.bytime.content[,6],type='l',lwd=3,col='black',ylab='',xlab='Republican News Exposure (Low to High)',las=1,ylim=c(-15,15),main='',cex.lab=1.2, cex.axis=1.2)
points(x=1:5,y=candidate.assessment.bytime.content[,6],pch=16,cex=1.5)
lines(x=1:5,y=candidate.assessment.bytime.content[,7],type='l',lwd=3,col='gray')
points(x=1:5,y=candidate.assessment.bytime.content[,7],pch=16,col='gray',cex=1.5)

dev.off()

###
#Regression Analysis
###

#Regress List of Survey Outcomes on a Treatment Varaible
	#Candidate Therm Difference
	#Candidate Choice
	#Immigration / Abortion / Refugees / Muslim / Sex Discrimination / Gov Spending / Insurance 
	#Out Party Marriage / Out Party Supporters Selfish / Out Party Supporters Narrow
	#Out Party Candidate Angry / Out Party Candidate Afraid / Out Party Candidate Disgusted

#Treatment Variables
	#Rating Content
	#Partisan Audience

treatment.variables <- c('rep.partisan','rating.average')
outcome.variables <- c('candidate.diff.change','candidate.trait.change','candidate.feel.change')
outcome.variables.post <- c('candidate.diff.post','candidate.trait.post','candidate.feel.post')
outcome.variables.pre <- c('candidate.diff.pre','candidate.trait.pre','candidate.feel.pre')
outcome.labels <- c('Therm','Trait','Emotion')

individual.scaled <- individual.use.election[,which(names(individual.use.election) %in% c(treatment.variables,outcome.variables,outcome.variables.pre,outcome.variables.post))]
individual.scaled <- as.data.frame(scale(individual.scaled, center = T, scale = T)) ## Standardize all non-binary covariates for easier comparison
individual.use.election <- cbind.data.frame(individual.use.election[,c('pid7_pre','educ_pre','race_pre','gender_pre','faminc_pre','ideo5_pre','age','political.interest.factor','total.visits','election.visits')],individual.scaled)

treat.outcomes <- as.data.frame(matrix(data=NA,ncol=6))
names(treat.outcomes) <- c('treatment','outcome','point','se','point.lagged','se.lagged')
	
for(j in 1:length(treatment.variables)){
	
	treatment <- treatment.variables[j]
	
	point.est <- NA
	se <- NA
	
	point.est.lagged <- NA
	se.lagged <- NA
	
	for(k in 1:length(outcome.variables)){
		
		outcome <- outcome.variables[k]
		post.outcome <- outcome.variables.post[k]
		pre.outcome <- outcome.variables.pre[k]

		out.form <- as.formula(paste(outcome,'~',treatment,'+ pid7_pre + gender_pre + race_pre + educ_pre + political.interest.factor + age + faminc_pre + ideo5_pre + total.visits'))
		out.mod <- lm(formula=out.form,data=individual.use.election)
		
		out.form.lagged <- as.formula(paste(post.outcome,'~',treatment,'+',pre.outcome,'+ pid7_pre + gender_pre + race_pre + educ_pre + political.interest.factor + age + faminc_pre + ideo5_pre + total.visits'))
		out.mod.lagged <- lm(formula=out.form.lagged,data=individual.use.election)
		
		point.est[k] <- out.mod$coefficients[2]
		se[k] <- coef(summary(out.mod))[2,2]
		
		point.est.lagged[k] <- out.mod.lagged$coefficients[2]
		se.lagged[k] <- coef(summary(out.mod.lagged))[2,2]
	}
	
	save.coef <- cbind.data.frame(rep(treatment,length(outcome.variables)),outcome.variables,point.est,se,point.est.lagged,se.lagged)
	names(save.coef) <- c('treatment','outcome','point','se','point.lagged','se.lagged')
	treat.outcomes <- rbind(treat.outcomes,save.coef)
}
treat.outcomes <- treat.outcomes[-1,]
treat.outcomes$upper <- treat.outcomes$point + 2*treat.outcomes$se
treat.outcomes$lower <- treat.outcomes$point - 2*treat.outcomes$se

treat.outcomes$upper.lagged <- treat.outcomes$point.lagged + 2*treat.outcomes$se.lagged
treat.outcomes$lower.lagged <- treat.outcomes$point.lagged - 2*treat.outcomes$se.lagged

##
#Plot the results
##

pdf(height=4,width=8,file='figure-e1.pdf')
par(mar=c(4,4.2,4,2),mfrow=c(1,2))
plot(y=length(outcome.variables):1+.1,x=treat.outcomes$point[1:length(outcome.variables)],pch=16,yaxt='n',ylab='',xlab='Republican Media Exposure',main="Site Audience",xlim=c(-.2,.2),ylim=c(.8,3.2),cex=1.4)
segments(length(outcome.variables):1+.1,x0=treat.outcomes$upper[1:length(outcome.variables)],x1=treat.outcomes$lower[1:length(outcome.variables)],lwd=3)
axis(side=2,at=length(outcome.variables):1,labels=outcome.labels,las=1)
abline(v=0,lty=2)
points(length(outcome.variables):1-.1,x=treat.outcomes$point.lagged[1:length(outcome.variables)],pch=16,col='gray',cex=1.4)
segments(length(outcome.variables):1-.1,x0=treat.outcomes$upper.lagged[1:length(outcome.variables)],x1=treat.outcomes$lower.lagged[1:length(outcome.variables)],lwd=3,col='gray')

plot(y=length(outcome.variables):1+.1,x=treat.outcomes$point[(length(outcome.variables)+1):(length(outcome.variables)*2)],pch=16,yaxt='n',ylab='',xlab='Republican Media Exposure',main="Article Slant",xlim=c(-.2,.2),cex=1.4,ylim=c(.8,3.2))
segments(length(outcome.variables):1+.1,x0=treat.outcomes$upper[(length(outcome.variables)+1):(length(outcome.variables)*2)],x1=treat.outcomes$lower[(length(outcome.variables)+1):(length(outcome.variables)*2)],lwd=3)
axis(side=2,at=length(outcome.variables):1,labels=outcome.labels,las=1)
abline(v=0,lty=2)
points(length(outcome.variables):1-.1,x=treat.outcomes$point.lagged[(length(outcome.variables)+1):(length(outcome.variables)+length(outcome.variables))],pch=16,col='gray',cex=1.4)
segments(length(outcome.variables):1-.1,x0=treat.outcomes$upper.lagged[(length(outcome.variables)+1):(length(outcome.variables)+length(outcome.variables))],x1=treat.outcomes$lower.lagged[(length(outcome.variables)+1):(length(outcome.variables)+length(outcome.variables))],lwd=3,col='gray')
dev.off()

sink()